///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file formic/wfn/coeff_obj.h
///
/// \brief   header file for the coefficient object class
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#ifndef FORMIC_COEFF_OBJ_HEADER
#define FORMIC_COEFF_OBJ_HEADER

#include<vector>
#include<string>

#include<boost/shared_ptr.hpp>

#include<formic/numeric/mantissa_exponent.h>

//#include<src/pcps.h>
//#include<src/input.h>
//#include<src/wfn_data.h>
//#include<src/configuration.h>

namespace formic {

  class Archive;
  class Configuration;
  class ConfigDiff;
  class InputBase;

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief  a class used to create simplifying typedefs for shared pointers
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template<class T> class PtrTypes {
    public:
      typedef boost::shared_ptr<T> p;
      typedef std::vector<boost::shared_ptr<T> > vp;
      typedef boost::shared_ptr<std::vector<boost::shared_ptr<T> > > pvp;
  };

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   Base class for an object representing a coefficient for a given occupation number
  ///          vector.
  ///
  /// The object should be capable of calculating the coefficient, local energy, overlap contraction,
  /// and Hamiltonian contraction for the current configuration as well as the ratio between the
  /// the current and a proposed configuration.
  ///
  /// The coefficient may depend on integer and floating point wavefunction data as well as other
  /// coefficient objects, allowing for the overall wavefunction coefficient to be a non-trivial
  /// combination of many different types of sub-coefficients (such as a Jastrow factor coefficient
  /// multiplied by a Slater determinant coefficient).
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template<class S> class CoeffObj {

    private:

      void not_implemented(const std::string & func) const;

    protected:

      /// \brief  the name of the coefficient
      std::string _name;

      /// \brief  the alpha indices of the configurations this process handles for exact sampling
      std::vector<int> _exact_ia_list;

      /// \brief  the beta indices of the configurations this process handles for exact sampling
      std::vector<int> _exact_ib_list;

//      /// \brief  the value of the coefficient for the currently initialized configuration
//      pcps::MantExp<S> _curr_val;
//
//      /// \brief  pointers to floating point wavefunction data on which the coefficient depends
//      typename pcps::PtrTypes<pcps::WfnData<S> >::pvp _wf_deps;
//
//      /// \brief  pointers to integer wavefunction data on which the coefficient depends
//      typename pcps::PtrTypes<pcps::WfnData<int> >::pvp _wi_deps;
//
//      /// \brief  pointers to other coefficient objects on which the coefficient depends
//      typename pcps::PtrTypes<pcps::CoeffObj<S> >::pvp _co_deps;

    protected:

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief updates the coefficient object's variables
      ///
      /// \param[in]      dfc      the derivative function coefficients.  size == this->n_der()
      /// \param[in]      alpha    the constant to scale the update by
      ///
      ///  |Psi(x_1, x_2, ...)> + sum_{i=1} dfc_i |dPsi/dx_i>
      ///     --->  |Psi(x_1 + alpha * dfc_1, x_2 + alpha * dfc_2, ...)>
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual void update_variables_detail(const S * const dfc, const S alpha) {
        this->not_implemented("update_variables_detail");
      }

      // constructors are protected as they are only meant for use by the childrens' constructors
      CoeffObj();
      CoeffObj(const CoeffObj<S> & other);

    public:

      /// \brief clears the lists of configuration indices for exact sampling, indicating that they will
      ///        will need to be re-computed before the next exact sampling is done
      void clear_exact_sample_vecs() {
        _exact_ia_list.clear();
        _exact_ib_list.clear();
      }

//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief constructor that initializes the pointers to point vectors of the specified size
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      CoeffObj(const int nwf = 0, const int nwi = 0, const int nco = 0)
//        : _wf_deps( new typename pcps::PtrTypes<pcps::WfnData<S>   >::vp(nwf) ),
//          _wi_deps( new typename pcps::PtrTypes<pcps::WfnData<int> >::vp(nwi) ),
//          _co_deps( new typename pcps::PtrTypes<pcps::CoeffObj<S>  >::vp(nco) )
//      {}

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief destructor
      ///
      /// Declared virtual so derived classes' destructors are called when deleting parent class's pointer.
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual ~CoeffObj() {}

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the name of the component
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      const std::string & name() const { return this->_name; }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the name of the child class
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual std::string class_name() const { this->not_implemented("class_name"); return std::string(); }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Initializes the object for the given lattice configuration.
      ///
      /// \param[in]      config  the lattice configuration we want the object initialized for
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual void initialize(const formic::Configuration & config) { this->not_implemented("initialize"); }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Updates the object so that the current configuration becomes the new configuration,
      ///        i.e. current = current + diff
      ///
      /// \param[in]      diff    an object describing the change to the configuration
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual void update(const formic::ConfigDiff & diff) { this->not_implemented("update"); }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the value of the coefficient for the current configuration
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual formic::MantExp<S> get_coeff() { this->not_implemented("get_coeff"); return formic::MantExp<S>(); }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the value of the coefficient for an arbitrary configuration
      ///
      /// \param[in]      config  the lattice configuration for which we want the coefficient
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual formic::MantExp<S> get_coeff(const formic::Configuration & config) { this->not_implemented("get_coeff(config)"); return formic::MantExp<S>(); }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the ratio of the coefficient at a new configuration (current + diff) to
      ///        that at the current configuration
      ///
      /// \param[in]      diff     the differece between the new and current configurations
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual S get_ratio(const formic::ConfigDiff & diff) { this->not_implemented("get_ratio"); return S(); }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the local energy for the current configuration
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual S local_energy() { this->not_implemented("local_energy"); return S(); }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the diagonal contribution to the local energy for the current configuration
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual S local_diag_energy() { this->not_implemented("local_diag_energy"); return S(); }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief gives the derivative ratios for the current configuration
      ///
      /// \param[out]     ratios     on exit, the derivative ratios.
      ///                            for ign_fixed == true, size == this->n_der()
      ///                            for ign_fixed == false, size == this->n_var()
      /// \param[in]      ign_fixed  whether to ignore variables that are fixed
      ///
      ///   The ratio for variable i is     <n|Psi_i> / <n|Psi>
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual void der_ratios(S * const ratios, const bool ign_fixed = true) { this->not_implemented("der_ratios"); }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the number of derivatives, which is the number of non-fixed variables
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual int n_der() const { this->not_implemented("n_der"); return 0; }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the number of variables (including "fixed" variables)
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual int n_var() const { this->not_implemented("n_var"); return 0; }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Returns the identity contraction with a derivative ratio vector.  Note that the first
      ///        element of the vector corresponds to the no-derivative term.
      ///
      /// \param[in,out]  dfc      the derivative function coefficients.  size == 1+this->n_der()
      ///
      /// \return  sum_i <n|I|Psi_i> dfc_i / <n|Psi>
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual S der_ic(const S * const dfc) { this->not_implemented("der_ic"); return S(); }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Returns the Hamiltonian contraction with a derivative ratio vector.  Note that the first
      ///        element of the vector corresponds to the no-derivative term.
      ///
      /// \param[in,out]  dfc      the derivative function coefficients.  size == 1+this->n_der()
      ///
      /// \return  sum_i <n|H|Psi_i> dfc_i / <n|Psi>
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual S der_hc(const S * const dfc) { this->not_implemented("der_hc"); return S(); }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Computes the local energy and its derivatives with respect to the wavefunction variables
      ///
      /// \param[out]     local_e  On exit, the local energy.
      /// \param[out]     le_ders  On exit, the derivative of the local energy.  size == this->n_der()
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual void local_energy_derivatives(S & local_e, S * const le_ders) {
        this->not_implemented("local_energy_derivatives");
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  Computes the orbital coeff matrix derivative ratios for the current configuration.
      ///         
      /// \param[in]      Cip      The orbital coefficient matrix.
      ///                          size == norb*norb
      /// \param[out]     dPsidCip On exit, the derivative ratios w.r.t. the orbital coefficient matrix.
      ///                          size == norb*norb
      ///
      ///         The orbital coefficient derivative ratio w.r.t. the (i,p) element of the orbital
      ///         coefficient matrix C is defined as  <dPsi/dC(i,p)|n> / <Psi|n>
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual void orb_coeff_der_ratios(const S * const Cip, S * const dPsidCip) {
        this->not_implemented("orb_coeff_der_ratios");
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  Gets the local configuration's weighted contribution to the energy derivatives w.r.t.
      ///         the one and two electron integral elements.
      ///         
      /// \param[in]      weight   Weight to give the local contributions.
      /// \param[out]     local_e  On exit, the local energy.
      /// \param[in,out]  mei_adj  Cumulative derivative adjoints for modified one electron integrals.
      ///                          size == norb*norb
      /// \param[in,out]  tei_adj  Cumulative derivative adjoints for two electron integrals.
      ///                          size == norb*norb*norb*norb
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual void local_energy_integral_der(const double weight, S & local_e, S * const mei_adj, S * const tei_adj) {
        this->not_implemented("local_energy_integral_der");
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief updates the coefficient object's variables
      ///
      /// \param[in]      dfc      the derivative function coefficients.  size == this->n_der()
      /// \param[in]      alpha    the constant to scale the update by
      ///
      ///  |Psi(x_1, x_2, ...)> + sum_{i=1} dfc_i |dPsi/dx_i>
      ///     --->  |Psi(x_1 + alpha * dfc_1, x_2 + alpha * dfc_2, ...)>
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void update_variables(const S * const dfc, const S alpha = formic::unity(S())) {
        this->update_variables_detail(dfc, alpha);
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief copies the coefficient's variables to the supplied vector
      ///
      /// \param[out]     vars     on exit, the variables. size == this->n_der()
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual void to_vector(S * const vars) { this->not_implemented("to_vector"); }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief copies the coefficient's variables from the supplied vector
      ///
      /// \param[in]      vars     the variables. size == this->n_der()
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual void from_vector(const S * const vars) { this->not_implemented("from_vector"); }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief reads in all the coefficient object's variables (even the "fixed" ones)
      ///
      /// \param[in]      vars     the variables. size == this->n_var()
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual void all_from_vector(const S * const vars) { this->not_implemented("from_vector"); }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns a vectorization of all the coefficient's variables (even the "fixed" ones)
      //
      /// \param[out]     vars     vector in which to store the variables. size == this->n_var()
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual void all_to_vector(S * const vars) const { this->not_implemented("all_to_vector"); }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief prints the coefficient's variables to the output file
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual void print_variables() { this->not_implemented("print_variables"); }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  updates the coefficient object with a rotation of the orbital basis
      ///
      /// \param[in]      norb     the number of orbitals
      /// \param[in]      orot     the orbital rotation matrix (new/old basis labelled by rows/cols).
      ///                          size == norb*norb
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual void apply_orbital_rotation(const int norb, const S * const orot) {
        this->not_implemented("apply_orbital_rotation");
      }

//      // ///////////////////////////////////////////////////////////////////////////////////////////////////
//      // /// \brief returns the value of the coefficient for the currently initialized configuration
//      // ///        with noise removed
//      // ///
//      // /// \param[in]      config   the current lattice and noise configuration
//      // ///
//      // ///////////////////////////////////////////////////////////////////////////////////////////////////
//      // virtual pcps::MantExp<S> noiseless_value(const pcps::Configuration<S> & config) { return this->_curr_val; }
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief returns the value of the coefficient for the currently initialized configuration
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      virtual pcps::MantExp<S> get_value() { return this->_curr_val; }
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief returns the value of the coefficient for an arbitrary configuration
//      ///
//      /// \param[in]      config   the configuration whose coefficient we want
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      virtual pcps::MantExp<S> get_value(const pcps::Configuration<S> & config) = 0;
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief returns the value of the coefficient for a configuration that differs from the
//      ///        currently initialized configuration as specified.
//      ///
//      /// \param[in]      new_config  the lattice configuration whose coefficient we want
//      /// \param[in]      old_config  the lattice configuration the object is currently initialized for
//      /// \param[in]      diff        an object specifying how the configurations differ
//      ///
//      /// \return the coefficient at the new configuration
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      virtual pcps::MantExp<S> get_value(const pcps::Configuration<S> & new_config,
//                                         const pcps::Configuration<S> & old_config,
//                                         const pcps::ConfigDiff<S> & diff) = 0;
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief Computes the ratios of the noiseless coefficients at each singly excited configuration
//      ///        to the noisy coefficient at the current configuration
//      ///
//      /// \param[in]      config      the lattice configuration the object is currently initialized for
//      /// \param[out]     a_ratios    on exit, the matrix of excitation ratios for alpha excitations
//      /// \param[out]     b_ratios    on exit, the matrix of excitation ratios for beta excitations
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      virtual void single_excitation_ratios(const pcps::Configuration<S> & config,
//                                            std::vector<S> & a_ratios,
//                                            std::vector<S> & b_ratios) = 0;
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief Adds the derivatives of the coefficient with respect to the wavefunction's data
//      ///        variables to the supplied vector after multiplying the derivatives by a constant.
//      ///        The derivatives are computed assuming that the object is currently initialized for the
//      ///        supplied lattice configuration.
//      ///
//      /// \param[in]      config  the lattice configuration the object is currently initialized for
//      /// \param[in]      alpha   constant to multiply the coefficient's derivatives by
//      /// \param[out]     der_vec vector to which to add the derivative contributions from this object
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      virtual void add_derivatives(const pcps::Configuration<S> & config,
//                                   const pcps::MantExp<S> & alpha,
//                                   std::vector<S> & der_vec) = 0;
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief For each orbital pair excitation, contracts the wavefunction's derivatives at the
//      ///        corresponding excited configuration with the supplied trial vector.
//      ///        It is assumed that the object is initialized for the supplied configuration.
//      ///
//      /// \param[in]      config      the lattice configuration to excite from.
//      /// \param[in]      alpha       constant to multiply the contraction by
//      /// \param[in]      trial_vec   the vector to be contracted
//      /// \param[in,out]  out_vec_a   vector to which to add the alpha orbital result
//      /// \param[in,out]  out_vec_b   vector to which to add the beta orbital result
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      virtual void contract_der_to_orb_pairs(const pcps::Configuration<S> & config,
//                                             const pcps::MantExp<S> & alpha,
//                                             const std::vector<S> & trial_vec,
//                                             std::vector<S> & out_vec_a,
//                                             std::vector<S> & out_vec_b) = 0;
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief For each wavefunction variable, contracts the various singly excited configurations'
//      ///        wavefunction derivatives with respect to that variable with the supplied trial vector.
//      ///        It is assumed that the object is initialized for the supplied configuration.
//      ///
//      /// \param[in]      config        the lattice configuration to excite from.
//      /// \param[in]      alpha         constant to multiply the contraction by
//      /// \param[in]      trial_vec_a   the vector to be contracted for alpha orbitals
//      /// \param[in]      trial_vec_b   the vector to be contracted for beta orbitals
//      /// \param[in,out]  out_vec       vector to which to add the result of the contraction
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      virtual void contract_orb_pairs_to_der(const pcps::Configuration<S> & config,
//                                             const pcps::MantExp<S> & alpha,
//                                             const std::vector<S> & trial_vec_a,
//                                             const std::vector<S> & trial_vec_b,
//                                             std::vector<S> & out_vec) = 0;
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief Multiply the supplied trial vector by the wavefunction derivative outer product matrix
//      ///        (the local contribution to the overlap matrix) for the given configuration.
//      ///        It is assumed that the object is initialized for the supplied configuration.
//      ///
//      /// \param[in]      config      the lattice configuration to excite from.
//      /// \param[in]      alpha       constant to multiply the contraction by
//      /// \param[in]      trial_vec   the vector to be contracted
//      /// \param[in,out]  out_vec     vector to which to add the result of the contraction
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      virtual void contract_der_to_der(const pcps::Configuration<S> & config,
//                                       const pcps::MantExp<S> & alpha,
//                                       const std::vector<S> & trial_vec,
//                                       std::vector<S> & out_vec) = 0;
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief Initializes the object for the given lattice configuration.
//      ///
//      /// \param[in]      config  the lattice configuration we want the object initialized for
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      virtual void initialize(const pcps::Configuration<S> & config) = 0;
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief Updates the object so that it is initialized for a new configuration that differs from
//      ///        the currently initialized configuration as specified.
//      ///
//      /// \param[in]      new_config  the lattice configuration we want to object to be initialized for
//      /// \param[in]      old_config  the lattice configuration the object is currently initialized for
//      /// \param[in]      diff        an object specifying how the configurations differ
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      virtual void update(const pcps::Configuration<S> & new_config,
//                          const pcps::Configuration<S> & old_config,
//                          const pcps::ConfigDiff<S> & diff) = 0;
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief returns the number of floating point variables that the component depends upon
//      ///
//      /// \param[in,out]  counted_dependencies  a set containing the names of wavefunction floating point
//      ///                                       data dependencies that have already been counted towards
//      ///                                       the total number of variables
//      ///
//      /// \return the number of floating point variables that the component depends upon
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      int nfvar(std::set<std::string> & counted_dependencies = std::set<std::string>()) const {
//
//        // zero the total
//        int retval = 0;
//
//        // for each of this object's wavefunction floating point data dependencies, add its number of
//        // variables if they have not been counted already
//        typename pcps::PtrTypes<pcps::WfnData<S> >::vp::const_iterator w;
//        for (w = this->_wf_deps->begin(); w != this->_wf_deps->end(); w++)
//          if ( counted_dependencies.count((**w).name()) == 0 ) {
//            retval += (**w).size();
//            counted_dependencies.insert((**w).name());
//          }
//
//        // add the variables from each coefficient object dependency
//        typename pcps::PtrTypes<pcps::CoeffObj<S> >::vp::const_iterator c;
//        for (c = this->_co_deps->begin(); c != this->_co_deps->end(); c++)
//          retval += (**c).nfvar(counted_dependencies);
//
//        // return the total
//        return retval;
//
//      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief returns the number of nonzero matrix elements in the sparse preconditioner
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual int n_pc_nonzero() { this->not_implemented("n_pc_nonzero"); }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief computes the index arrays for the sparse preconditioner
      ///
      /// \param[in,out]  pc_n_inds  size 1+this->n_der().
      ///                            On exit, the number of nonzero matrix elements for each row.
      /// \param[in,out]  pc_inds    size this->n_pc_nonzero().
      ///                            On exit, the column indices of each rows' nonzero elements.
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual void set_pc_inds(int * const pc_n_inds, int * const pc_inds) {
        this->not_implemented("set_pc_inds");
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief clones the object
      ///
      /// \return  a shared pointer to a new clone of the object
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual boost::shared_ptr<formic::CoeffObj<S> > clone() const { this->not_implemented("clone"); }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief converts the object into the command that will create the object
      ///
      /// \param[in,out]  arch     an archive in which to store the command
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual void to_command(formic::Archive & arch) const { this->not_implemented("to_command"); }

      void prep_for_exact_sample(const formic::InputBase & userinp);

      /// \brief  return the alpha indices of the configurations this process handles for exact sampling
      const std::vector<int> & exact_ia_list() const { return _exact_ia_list; }

      /// \brief  return the beta indices of the configurations this process handles for exact sampling
      const std::vector<int> & exact_ib_list() const { return _exact_ib_list; }

  };

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  overload the << operator to enable easy archiving of coefficient objects
///
/// \param[in,out]  arch     an archive in which to store the coefficient object
/// \param[in,out]  c        the coefficient object we wish to store
///
/// \return the archive
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class T>
inline formic::Archive & operator<<(formic::Archive & arch, const formic::CoeffObj<T> & c) { c.to_command(arch); return arch; }

#endif
